// display the convolved momentum distribution
// the input file format is mom str tot
// by Asia Sun, Jun. 23, 2022
// Last modified: Apr. 2, 2023

static int co = 0; // graph color
void primp(TGraph *g); // primp the graph
const string dir = "";
// C15-n c.s. on C12 target @237MeV/u
const double csExp = 107.7019, dcsExp = 11.3779; // -n c.s.
const double csExp2n = 59.6636, dcsExp2n = 7.2530; // -2n c.s.

void momConv(){
  string file = dir+"mom.txt";
  TGraph *gtot = new TGraph(file.c_str(), "%lg %*lg %lg"); primp(gtot);
  TLegend *lg = new TLegend(0.6, 0.65, 0.89, 0.85); lg->SetBorderSize(0);
  lg->AddEntry(gtot, "Total");

  // the exp. momentum distribution //
  TFile *f = new TFile("ppl.root");
  TH1F *hC14 = (TH1F *)f->Get("hC14NET"); hC14->SetLineColor(4);
  TH1F *hC13 = (TH1F *)f->Get("hC13NET");
  hC14->Scale(csExp/hC14->Integral("width")); // scale to the total c.s.
  hC13->Scale(csExp2n/hC13->Integral("width")); // scale to the total c.s.
  lg->AddEntry(hC14, "^{15}C-n Mom. distri.");
  lg->AddEntry(hC13, "^{15}C-2n Mom. distri.");

  // scaling //
  // energy shift
  const double dm = gtot->GetMean() - hC14->GetMean()-10;
  cout << "gtot->GetMean(): " << gtot->GetMean() << endl;
  cout << "hC14->GetMean(): " << hC14->GetMean() << endl;
  cout << "dm: " << dm << endl;
  // amplitude scaling
  const double cce = hC14->Integral("width")/gtot->Integral();
  // const double cce = 1.;
  cout << "gtot->Integral(): " << gtot->Integral() << endl;
  cout << "hC14->Integral(): " << hC14->Integral() << endl;
  cout << "cce: " << cce << endl;
  const int n = gtot->GetN(); double x, y;
  for(int i = 0; i < n; i++){
    gtot->GetPoint(i, x, y); gtot->SetPoint(i, x-dm, y*cce);
  } // end for over i

  TCanvas *c = new TCanvas("c1", "c1", 800, 600);
  gtot->SetLineWidth(2); gtot->SetLineStyle(1);
  gtot->SetTitle("");
  gStyle->SetOptStat(0);
  TAxis *ax = hC14->GetXaxis(); ax->SetTitle("p_{c} [MeV/c]"); ax->CenterTitle();
  TAxis *ay = hC14->GetYaxis(); ay->SetTitle("d\\sigma_{-n}/dp_{c} [mb/(MeV/c)]");
  ay->CenterTitle(); ay->SetTitleOffset(1.30);
  hC14->SetTitle("");
  hC14->Draw();

  gtot->Draw("same");
  gStyle->SetOptStat(0);
  hC13->Draw("same");
  hC13->Fit("gaus","R", "", 8000, 9300);

  lg->Draw();
  c->Print("mom.pdf");
  c->Print("mom.root");
} // end main function mom

// primp the graph
void primp(TGraph *g){
  TAxis *ax = g->GetXaxis(); ax->SetTitle("p_{c} [MeV/c]"); ax->CenterTitle();
  TAxis *ay = g->GetYaxis(); ay->SetTitle("d\\sigma_{-n}/dp_{c} [mb/(MeV/c)]");
  ay->CenterTitle(); ay->SetTitleOffset(1.30);
  co++;
  if(5 == co || 10 == co) co++;
  g->SetLineColor(co); g->SetMarkerColor(co);
} // end function primp
